% GAUTIER LE BIHAN - 2020
% Replication files for "Shocks vs Menu Costs: Patterns of Price Rigidity in an Estimated Multi-Sector Menu-Cost Model?" Review of Economics and Statistics
%
% Table 6
clear;
tic

addpath('..\..\Utilities')  
load actual_moments_k



%ENTER HERE THE MAX NB OF SECTOR
max_sectors=10;
load actual_moments_k
weight_j=actual_moments_k(:,12)/100;

load ..\..\Estimation_param\MS_produits_CPlus\param



param=param;
%*******************************************STAT SECTEURS PARAM
coicop1=actual_moments_k(:,1);
weight=actual_moments_k(:,12);
sum(weight);
sample_all=[coicop1 weight param];

moments_mean=(sample_all(:,9:end)'*weight)./sum(weight);
size(sample_all)

%/************* MEDIAN ALL PARAMETERS
%%%%%%%%%%%%%%%%%%%%%%TOTAL
p0=abs(tanh(sample_all(:,3)));
mu_c=abs(sample_all(:,4));
sig_a=exp(sample_all(:,5));
rho_a=ones(227,1)*0.6946;

param_all=[p0 mu_c sig_a rho_a];
size(param_all)
param_all=param_all(mu_c<0.4,:);
size(param_all)

weight=sample_all(mu_c<0.4,2);
v=size(param_all,1);
weight_all=weight;

wmean_all=(param_all'*weight_all)./sum(weight_all)
%wmean_all(2)=0.04*4/3

for i=1:size(param_all,2)
        [sortx,order] = sort(param_all(:,i));
        sortw = weight_all(order);

         midpoint = sum(sortw)/2;
         csumw = cumsum(sortw);
         j = find(csumw<=midpoint,1,'last');
         j
            if csumw(j) == midpoint
              wmedwall(:,i) = mean(sortx([j j+1]));
            else
                wmedwall(:,i) = sortx(j+1);
            end
         
        end  

 



size_p=size(param_all,1);
param=param_all;
for jj=2:max_sectors;
    number=jj;
%figure;hist(param_sect(:,3));

percentiles=ones(number-1,5);
    for i=1:size(param,2)
        [sortx,order] = sort(param(:,i));
        sortw = weight(order);
        
        for n=1:number-1
         midpoint = n*sum(sortw)/number;
         csumw = cumsum(sortw);
         j = find(csumw<=midpoint,1,'last');
         
        if isempty(j)==1
             j=1;
         end
            if csumw(j) == midpoint
              percentiles(n,i) = mean(sortx([j j+1]));
            else
                
                percentiles(n,i) = sortx(j+1);
            end
     percentiles(n,end)= csumw(end);
        end
    end  
    all_per(jj-1).all_per=percentiles;
end 
IRFd(5).IRFd={};



for ii=1:3

    stat_aggall=zeros(8,7)
    stat_sect=zeros(8,30)
    IRF_all=zeros(201,8)
for zz=1:8
    

%freq=moments(:,1);
choix_param=ii;
sig=param(:,choix_param);
quantile=zeros(size_p,max_sectors-1);

for jj=1:max_sectors-1

per_sig=all_per(jj).all_per;
per_sig=per_sig(:,choix_param);
var_dum=zeros(size_p,jj);
    for n=1:jj
    var_dum(:,n)=(sig>=per_sig(n));
   
    end
    if jj==1
    quantile(:,jj)=var_dum;
    else
     quantile(:,jj)=sum(var_dum')';
    end
end


for kk=1:max_sectors-1

quantile_k=quantile(:,kk);
max(quantile_k)
     wmean_quant=zeros(max(quantile_k)+1,size(param,2));
     weight_quant=zeros(max(quantile_k)+1,1);
    wmedw=zeros(max(quantile_k)+1,size(param,2));
    for k=1:max(quantile_k)+1
    
     param_sect=param((quantile_k==k-1), :);
     v=size(param_sect,1);
     weight_sect=weight((quantile_k==k-1),:);
    
     weight_quant(k,:)= sum(weight_sect);
     wmean_quant(k,:)= (param_sect'*weight_sect)/weight_quant(k,:);
     
     %***** MEDIAN by quintile****
     
       for i=1:size(param_sect,2)
        [sortx,order] = sort(param_sect(:,i));
        sortw = weight_sect(order);

         midpoint = sum(sortw)/2;
         csumw = cumsum(sortw);
         j = find(csumw<=midpoint,1,'last');
         %j
            if csumw(j) == midpoint
              wmedw(k,i) = mean(sortx([j j+1]));
            else
                wmedw(k,i) = sortx(j+1);
            end
         
     end  

end


 param_ms(kk+1).param_ms=wmean_quant;
  weight_ms(kk+1).weight_ms=weight_quant;
  param_msmed(kk+1).param_msmed=wmedw;
end;

 param_ms(1).param_ms=[wmean_all'];
  weight_ms(1).weight_ms=57.0181;
  param_msmed(1).param_msmed=[wmedwall]
weight_ms
s=40;
%stat_mean=zeros(s,7);
%for jj=[1 2 3 4 5 6 7 9 11 13 15 17 19];
 % for jj=[3 4 5 6 7 9 11 13 15 17 19];
    tt=5;
   
    param_med=param_ms(tt).param_ms;
    param_med=param_med((param_med(:,1)~=0),:)
    param_med
   
   if zz==1
   param_med(:,1)=ones(1,tt)*wmean_all(1);
   %param_med(:,3)=ones(1,tt)*wmean_all(3);
   end
    
   if zz==2
   param_med(:,2)=ones(1,tt)*wmean_all(2);
   %param_med(:,3)=ones(1,tt)*wmean_all(3);
   end

   if zz==3
   param_med(:,3)=ones(1,tt)*wmean_all(3);
  % param_med(:,2)=ones(1,tt)*wmean_all(2);
   end
   
   if zz==4
   
   end
    
   if zz==5
   param_med(:,1)=ones(1,tt)*wmean_all(1);
   param_med(:,2)=ones(1,tt)*wmean_all(2);
   param_med(:,3)=ones(1,tt)*wmean_all(3);
   end

   if zz==6
   param_med(:,2)=ones(1,tt)*wmean_all(2);
   param_med(:,3)=ones(1,tt)*wmean_all(3);
   end
    
   if zz==7
   param_med(:,1)=ones(1,tt)*wmean_all(1);
   param_med(:,3)=ones(1,tt)*wmean_all(3);
   end
   
   if zz==8
   param_med(:,1)=ones(1,tt)*wmean_all(1);
   param_med(:,2)=ones(1,tt)*wmean_all(2);
   end
param_med
   %param_med(ii,:)=ones(4,1)*param_bar(ii);
    
    weight_sect=weight_ms(tt).weight_ms;
    weight_sect=weight_sect((weight_sect(:,1)~=0),:)
    weight_sect
    p0=param_med(:,1)';
    mu_c=param_med(:,2)';
    sig_eps_a=param_med(:,3)';
    rho_a=ones(1,tt)*0.6946;   

     w=[weight_sect'];
    sm_m=0;
    param_e=[p0 mu_c sig_eps_a rho_a];

   [avmoments, moments, VC, IRF]=geNCalvoPlus_2Agg(sm_m, p0, mu_c, sig_eps_a, rho_a, w);
    NSect = size(w,2);


    stat_sect(zz,1:tt*6)=[ moments];
   % save stat_sect stat_sect;

  
    

end
  

      
  st_sector(ii).st_sector= stat_sect;
  save st_sector st_sector;
   
  

  

  
end 
  